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Experimental data of the DNA cyclization (J-factor) at short length scales, as a way to study 
the elastic behavior of tightly bent DNA, exceed the theoretical expectation based on the wormlike 
chain (WLC) model by several orders of magnitude. Here, we propose that asymmetric bending 
rigidity of the double helix in the groove direction can be responsible for extreme bendability of 
DNA at short length scales and it also facilitates DNA loop formation at these lengths. To account 
for the bending asymmetry, we consider the asymmetric elastic rod (AER) model which has been 
introduced and parametrized in an earlier study [Eslami-Mossallam, B.; Ejtehadi, M. R. Phys. Rev. 

E 2009, 80, 011919]. Exploiting a coarse grained representation of DNA molecule at base pair (bp) 
level, and using the Monte Carlo simulation method in combination with the umbrella sampling 
technique, we calculate the loop formation probability of DNA in the AER model. We show that, 
for DNA molecule has a larger J-factor compared to the WLC model which is in excellent agreement 
with recent experimental data. 


I. INTRODUCTION 

Studying the elastic behavior of DNA molecules is im¬ 
portant for understanding its biological functions. One of 
the most popular theoretical models to explain the elas¬ 
tic behavior of DNA is the harmonic elastic rod model, 
also called the wormlike chain (WLC) model [I] [2]. In 
this model it is assumed that the elastic energy is a har¬ 
monic function of local deformations. The WLC model 
can predict very accurately the elastic properties of long 
DNA molecules and yielding a persistence length of about 
50 nrn for DNA [2 [3| . However, recent experimental data 
suggest that, short DNA molecules are much more flexi¬ 
ble than what is predicted by the WLC model [2HE1- For 
example, loop formation probability, i.e. the J-factor [3], 
for DNA molecules shorter than 100 bp (~ 34 nm) is sev¬ 
eral orders of magnitude higher than the prediction of the 
WLC model [21 Ei sSJ HD]. 

Different experimental procedures have been used 
to measure the cyclization probability for short DNA 
molecules. For example, in Cloutier and Widom’s 
work J5] the DNA molecules have short sticky ends. 
Therefore, when the two DNA ends are close to each 
other torsional and axial alignment are required to form 
a DNA loop, which is then stabilized by the ligase. Thus 
the J-factor depends on the concentration of the ligase in 
the experiment m On the other hand, Vafabakhsh and 
Ha have used DNA molecules with long sticky ends [8]. 
In this case it is expected that the rate of loop forma¬ 
tion depends only on the probability that the two DNA 
ends reach together, and thus is directly related to DNA 
elasticity. In the both experiments the persistence length 
of short DNA molecules appear to be much shorter then 
50 nm. 

It has been suggested that the anomalous elastic be¬ 
havior of short DNA molecules is a consequence of forma- 
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tion low-energy kinks in highly bent DNA molecules jTTE 
uni. Also it has shown that local DNA melting of 
the cyclized DNA increases the J-factor at short length 
scales m S3- Also, there has been efforts to ex¬ 
plain the effect by introducing more structural details 
to elastic model (e.g., considering cooperativty mm 
or anisotropy mm in bending rigidity of DNA). But 
these efforts were not successful, as it has been shown 
that anisotropy has no significant effect in these dimen¬ 
sions [2U] and even leads to the stiffening of DNA if the 
molecule is confined in a two dimensional surface ED- 
The effect of electrostatic interaction has been studied 
which can increase the j-factor [22], but in this study we 
assume the DNA molecule is neutral and the electrostatic 
interaction would be ignored. 


The DNA molecule in its B form suggests that DNA 
has an asymmetric structure, in the sense that the op¬ 
posite grooves of the DNA are not equal in size and the 
structure [23]. Thus, one expects that the energy re¬ 
quired to bend the DNA is not only anisotropic, but also 
asymmetric, i.e. the energy required to bend the DNA 
over its major groove is not equal to the energy needed 
to bend over its minor groove to the same angle. There 
are theoretical analysis [22], experimental evidences [25] . 
and simulation studies [26] which show that the bending 
asymmetry may affect the elastic behaviour of DNA. In 
previous work we have introduced the asymmetric elastic 
rod (AER) model to account for the asymmetric bend¬ 
ing of DNA 127. . In this paper, we evaluate the elastic 
properties of the AER model, namely the DNA loop¬ 
ing probability, to reveal the relevance of the asymmetric 
bending for short extremely bent DNA molecules. To 
this end we exploit the Monte Carlo (MC) simulation in 
combination with the umbrella sampling (US) technique, 
which enables us to efficiently sample the rare cyclization 
events. We show that the AER model is in an excellent 
agreement with the experimental data of the J-factor at 
short length scales. 
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II. MODEL AND METHOD 
A. Asymmetric elastic rod model 

In the elastic rod model, a DNA is represented by a 
flexible inextensible rod mm, which can be deformed in 
response of the external forces or torques. Here we use 
the discrete elastic rod model El Ell, where the rod is 
discretized into segments each representing a DNA base 
pair. In this model, the internal degrees of freedom of 
the base pairs are neglected, and each base pair is consid¬ 
ered as a rigid body. A local coordinate system (material 
frame) with an orthonormal basis {e?i, d 2 , d 3 } is attached 
to each base pair. As depicted in Figure [TJ d 3 is perpen¬ 
dicular to the base pair surface, d\ lies in the base pair 
plain and points toward the major groove, and d 2 is de¬ 
fined as d 2 = d 3 x di. Since it is assumed that the DNA 
is inextensible, each base pair only has three rotational 
degrees of freedom, and the position of the (fc + l)th base 
pair with respect to the fcth base pair is denoted by the 
vector r(fc) which is given by Ei 

r(k) = l 0 d 3 (k) , (1) 

where Iq = 0.34 nrn is the base pair separation. The ori¬ 
entation of the (fc + l)th base pair with respect to the 
fcth base pair is determined by a rotation transforma¬ 
tion R(fc), which can be parametrized by a vector 0(fc). 
The direction of 0(fc) is normal to the plain of rotation 
of the fcth base pair, and its magnitude determines the 
rotation angle. The components of 0(fc) in the local coor¬ 
dinate system attached to the fcth base pair are denoted 
by 0i(fc), 02 (fc), and © 3 (fc). These components can be 
regarded as three rotational degrees of freedom of the 
base pairs around c?i, <i 2 and d 3l and are called tilt, roll, 
and twist respectively [25] , If the values of these three 
angles are known for all base pairs, the conformation of 
the DNA can be uniquely determined. 



FIG. 1: Parametrization of the elastic rod. The local frame 
{di, d 2 , d 3 } is attached to the rod. 

For an inextensible DNA with N base pair steps, the 
elastic energy depends on the spatial angular velocity 


f2(fc) = pp, then the elastic energy of the AER model 
m can be written as 


E = 


N 




AS 

k 


[H(fc)], 


( 2 ) 


with 

4 as Pm = 

k B T 


-Ai fli(fc)" + —A 2 ^2(fc ) 2 + 

-C (0 3 (fc) — w 0 ) 2 + F 2 fl 2 (fc) 3 + 


-G 3 (0 1 (fc) 4 + 0 2 (fc) 4 ) 


4! 


(3) 


where k B is the Boltzmann constant and T is the ab¬ 
solute temperature. The first three terms in equation 
([ 3 ]) correspond to the harmonic part of the elastic en¬ 
ergy, which also appear in the WLC model. Ai, A 2 and 
G are the harmonic elastic constants of DNA, and Uq 
is its intrinsic twist. The remaining terms in equation 



FIG. 2: (Color online) £2 as a function of for two different 
models, “W” and “A” (Table |T|. The symmetric model “W”, 
dashed (black) curve, has one minimum and its curvature is 
positive everywhere. The asymmetric model “A”, solid (red) 
curve, has two minima. £\ as a function of fii for both the 
models remain always convex (inset). 

([3]) constitute the anharmonic parts of the elastic en¬ 
ergy. The term +1/3! F 2 fl 3 accounts for the asymmetric 
structure of DNA in the bending energy, while the term 
jt G 3 (fli(fc) 4 + fl 2 (fc) 4 ) preserves the stability and the 
consistency of the model. 

Since there is no coupling term in the model, roll, tilt, 
and twist can be regarded as independent deformations, 
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and the energy density can be decomposed into three 
separate terms 

£ AS [fir, n 2 , n 3 ] = fi AS [^i] + £ 2 AS p 2 ] + £ AS [^ 3 ], (4) 


where 


£ AS [fil] = k B T n? + ^ g 3 Qi 


ASrn_i _ u-rr\l M ^2 _L ^2 ^3 + J_ g 3 ^4 


^2 " [^ 2 ] = feT |^ri 2 Si 2 ' g | 


4! 


and 

£ 3 as [^ 3 ] = \k B rc{n 3 - Wo ) 2 . 


(5) 

( 6 ) 

(7) 


Here, we use parameters of the model E3, which is ob¬ 
tained by fitting the AER model to the experimental 
data of Wiggins et al. [6] (see the first row of Table [T] 
as model “A”). In this parametrization, the bending 
anisotropy also is considered, where Ai 7^ A 2 . With 
this parametrization the roll energy, £ 2 , has two minima 
(solid, red curve in Figure [2]); one at fl 2 = 0 and the 
second one Sl 2 = — 3.3 nm _T corresponds to a negative 
roll of about 60°, with a roll energy about 8 k B T. The 
energy barrier between two minima is about 9 k B T. The 
existence of a second minimum in £ 2 can lead to the for¬ 
mation of kinks in the minor groove direction. With a 
large energy barrier between the two minima, one ex¬ 
pects that the kinks rarely form in a free DNA at room 
temperature. However, if the DNA is forced to adopt a 
tightly bent conformation the probability of kink forma¬ 
tion increases significantly. 

The possibility of kink formation in the DNA struc¬ 
ture has been considered previously by other authors. 
An atomistic structure for a kinked DNA has been pro¬ 
posed by Crick and Klug [24] . who suggest that DNA 
can form a kink in the minor groove direction. Also, 
molecular dynamics simulations on a 94 bp minicircle |26l 
show that kinks are formed, with the same structure pre¬ 
dicted by Crick and Klug. A simple model has been 
presented by Nelson, Wiggins, and Phillips to describe 
the elasticity of kinkable elastic rods jT2] . This model is 
mathematically equivalent to the models of local DNA 
melting [141 FTSl 29]. Recently, Vologodskii and Frank- 
Kamenetskii have proposed another model for the kink 
formation in DNA m- In all of these models, the kinks 
are isotropic, i.e. they can be formed in any direction 
with equal probability. On the contrary, in the AER 
model there is a privileged direction for the kink forma¬ 
tion, i.e. the groove direction. 

In order to compare the AER model with the WLC 
model, we also use another set of parameters here, which 
are given in second row of Table [T] as model “W”. As we 
will show in result section, at long length scales, these two 
models are equivalent and they yield the same persistence 
length, l p = 51 nm. The tilt and roll energies, £\ and £ 2 , 
of these two models are compared with each other in 
Figure [2] 


model 

Ai (nm) 

A 2 (nm) 

F (nm) 

G (nm) 

C (nm) 

cu 0 (1/nm) 

A 

87.00 

43.50 

7.90 

3.20 

100 

1.8 

W 

51.00 

51.00 

0.0 

0.0 

100 

1.8 


TABLE I: Two different parameter sets for the AER and 
WLC models, indicated as models “A” and “W” respectively 
in this study. The parameter set “A” is obtained from the 
fitting of the model to the experimental data of Wiggins et 
al. m and its effective persitence length is about 51 nm, 
same as molde “W”. 


B. Calculation of J-factor 


The loop formation probability of a polymer, which is 
known as J-factor GEO. is defined as the probability that 
the two ends of the polymer meet each other with axial 
and torsional alignment. For simplicity we neglect the 
torsional alignment and only require that the two ends 
are close to each other while the two terminal tangent 
vectors are parallel. Denoting the separation between the 
two ends by r, and the angle between the two terminal 
tangent vectors by 9 = cos -1 (d 3 (0) • d 3 {L )) (see figure [3J , 
the J-factor for a DNA of length L in molar unit is given 

by [501131! 



FIG. 3: (Color online) Schematic figure of a chain with termi¬ 
nal tangent vectors, d 3 (0) and d 3 (L), which are indicated by 
red and blue arrows, respectively. The end-to-end distance is 
shown by r, and 9 is the angle between the tangent vectors. 
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J(L)= lim —— -—., 
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K{r, L) dr 
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where 7 = cosd, N a is Avogadro’s Number, K(r,L) is 
the normalized distribution function of r, and P 7 (7, L) is 
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the normalized distribution function of 7 under condition 
of r = 0 . 

The bending energy of DNA depends on the bending 
direction. Thus there is a implicit bend-twist coupling in 
this model. However this coupling can barely affect the 
J-factor if the DNA length is much larger than the heli¬ 
cal pitch. Therefore we expect that for both models, the 
average trend of the J-factor is given by the equation ([ 8 ]). 
For end distances near zero ( r/L m 0), the radial distri¬ 
bution function, A\'(r, A), is proportional to r 2 

A'(r, A)— >4:Trr 2 K 0 (L ) for r/L —> 0, (9) 

where Kq(L) is a r-independent function [32] , In addition 
we have 

1 r 1 

iim -- / P 7 ( 7 , L) dy = P 7 (7 = 1, A), (10) 

70^1 I -70 J l0 

then the equation ([ 8 ]) can be written as 

J(A) = J 0 (A)x2P 7 (1,A) (11) 

where 

J 0 (L) = —pA'o(A), (12) 

1 'a 

is the unconstrained J-factor which dose not involve axial 
and torsional alignment. 

C. Simulations 

We exploited a Metropolis Monte Carlo (MC) simu¬ 
lation to evaluate the statistical properties of DNA. We 
do not include the self avoiding in the simulations, since 
the probability of self crossing is small for the short simu¬ 
lated DNA molecules. For short DNA molecules the loop 
formation probability is very low, and thus the DNA cy- 
clization events are too rare to be observed in the simu¬ 
lations. To overcome this problem, we used the method 
of Umbrella sampling (US) [33] to evaluate the distribu¬ 
tion functions K{r, A) and P 7 ( 7 , A). To calculate K(r , A) 
the reaction coordinate is the end-to-end distance, r, we 
divided r into 100 successive windows, and for each win¬ 
dow performed a separate MC simulation in which a har¬ 
monic bias potential is applied to the end-to-end distance 
of the DNA. All the harmonic potential have a common 
spring constant A/ = 0.3fcBA/nm 2 , and the minimum 
of each potential lies at the center of the corresponding 
window. We then found the biased distribution for each 
individual simulation and used the Weighted Histogram 
Analysis Method (WHAM) [34] to reconstruct the unbi¬ 
ased distribution function. To calculate P 7 ( 7 , A), we set 
the end-to-end distance to zero, then perform another 
US by dividing the range of variation of 7 into 100 win¬ 
dows and applying a harmonic bias potential with spring 
constant of A * = 40 k B T, in each window. The unbiased 
distribution function P 7 ( 7 , A) is then found by WHAM. 


In each individual simulation during the umbrella sam¬ 
pling procedure, the first 10 5 MC steps were disregarded 
to ensure the equilibration of the system and the next 
2 x 10 6 MC steps were considered for sampling. We per¬ 
form 5 independent simulations in each window to esti¬ 
mate errorbars. 



FIG. 4: (Color online) Monte Carlo simulation results of the 
effective bending energy, E e g, as a function of the bending 
angle, 9, for various chain lengths A = 3.4 11 m (black), A = 
6 .8nm (red), A = 13.6nm (blue), A = 27.2nm (green) and 
A = 51.0 nm (magenta). Solid and dashed curves correspond 
to “A” and “W” models, respectively (see Table [I]. 


III. RESULTS AND DISCUSSION 

A. Persistence Length and Effective Bending 
Energy 

For a long DNA of length A, the persistence length, l p , 
is defined as [131 . 

(cos(0)) = exp(— L/l p ). (13) 

With the parameters given in Table [I] the asymmetric 
model (model “A”) and the wormlike chain model (model 
“W”) have a common persistence length of about 51 nm. 
This means that for long DNA molecules with small de¬ 
formations the two model are equivalent, and thus the 
asymmetric model is effectively reduced to an isotropic 
wormlike chain model. On the other hand, at large bend¬ 
ing angles, one expects that the asymmetric structure of 
DNA affects its elastic properties. To show this, we eval¬ 
uated the effective bending energy as a function of the 
bending angle 9, which is defined as [351 

Al eff (0,A) = -fc B Tln^^, (14) 

suit/ 
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where P(0, L ) is the distribution function of the bending 
angle 9 of a DNA of length L. Figure [4] compares the ef¬ 
fective energies of the AER model (solid curves) with the 
WLC model (dashed curves) for different DNA lengths. 
One can see that, at small bending angles, both models 
follow a common parabola. However, at large bending 
angles, the effective bending energy of the asymmetric 
model falls beneath the parabola, which leads to extreme 
bendability of DNA or formation of kinks J22- The ef¬ 
fect is suppressed as the DNA length increases. It is well 
expected that, for long enough DNA, the effective energy 
is independent on the structural details and it converges 
to a parabola [5] . 

As Figure [4] shows, the transition between the har¬ 
monic and non-harmonic region is smooth. This is be¬ 
cause in the AER model DNA preserves its resistance 
against bending even in kink conformation. In simpler 
versions of kinkable elastic rod models [H [13] , where 
the kinks are assumed to be completely flexible, there is 
a sharp transition in the curve of the effective bending 
energy between a parabola and a straight line with zero 
slope )I2]. 



0 10 20 30 

r (nm) 


FIG. 5: (Color online) The radial distribution (a) and free 
energy (b) as a function of end-to-end distance, r, for a DNA 
of length L = 36.66 nm (= 99 bp). Data points represent MC 
simulation results, where triangles (blue) and squares (red) 
correspond to the models “A” and “W”, respectively. Solid 
(black) curve corresponds to the theoretical prediction of the 
WLC model with l p = 51 nm [36] . Error bars (not shown) are 
about the size of the symboles. 


B. The end-to-end distribution functions 

Figures |5j(a) and |5](b) compare the radial distribution 
function, K(r,L) for a DNA of length L = 33.66 nm(= 
99bp). The triangles (blue) and squares (red) show MC 
simulation results for the models “A” and “W”, respec¬ 
tively (Table [I]). The solid (black) curve corresponds to 
the theoretical treatment of Samuel and Sinha [36] for 
the WLC model, which perfectly matches the simulation 
data. As can be seen in Figure [5ja) there is no significant 
difference between the two models at large end-to-end 
distances, while at short end-to-end distances the radial 
distribution function in the AER model significantly de¬ 
viates from that of the WLC model. 



r/L 



FIG. 6: (Color online) a) K(r, L) / Anr 2 as a function of r 
for different DNA lengths. Solid and dashed curves corre- 
sponde to the models “A” and “W” in Table [I] respectively, 
b) Kq / K™ as a function of DNA length. 

Figure [5](b) shows — ln( K(r, L )), the free energy of the 
DNA as a function of its end-to-end distance. The po¬ 
sition of the minimum in the free energy curve corre¬ 
sponds to the most probable end-to-end separation. A 
relaxed DNA molecule which is shorter than the persis¬ 
tence length tends to be almost straight. As can be seen 
in Figure [5](b), for L = 99 bp the free energy minimum 
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is very close to the total length of DNA. In this case we 
found that the average and the variance of the end-to-end 
distance in the WLC and AER models differ by less than 
0.1 percent. Therefore in the experiments which involve 
long free DNA molecules, such as DNA stretching experi¬ 
ment, the two models are indistinguishable. On the other 
hand, the asymmetric bending can significantly affect the 
outcome of the experiments performed on short, tightly 
bent DNA molecules, such as the DNA cyclization. 


Figure 6(a) compares the distribution functions, 
K(r,L)/4nr' z , of the model “A” (solid curves) with the 
model “W” (dashed curve) for different lengths L = 
16.66, 36.66, 50.66, 67.66, 84.66, 101.66 nm (49, 99, 149, 
199, 249 and 299 bp, respectively). One can see that 
the difference between the two models is disappeared as 
the DNA length increases. As expected for small end- 
to-end distances A'(r, L)/47rr 2 converges to a constant 
K 0 {L) (see equation which is proportional to the 
unconstrained J-factor Jo (A) (equation ©)• To calcu¬ 
late Kq(L) we average K(r, L) / Airr 2 in vicinity of r = 0 
over the range of 0 < r < ro, whe re rp is chosen such 


6(b) 


shows the ratio 


ag K(r 0 ,L)-4*r 0 Ko(L) _ Q Q1 Figure 

K£/I<w as a function of DNA length, were the super¬ 
scripts “A” and “W” refer to the AER and WLC models 
as parametrized in Table [ij As can be seen, while K A 
is several order of magnitude larger than I\^ for short 
DNA molecules, the ratio K A /K^ approaches unity as 
the DNA length increases. 
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FIG. 7: (Color online) Top: different schematic configura¬ 
tions of the chain for different angles between the two DNA 
tangents. Bottom: the MC simulation results for P 7 (7, L) as 
a function of 7 for a DNA loop of length 36.66nm(= 99bp). 
Triangles (blue) and squares (red) correspond to the models 
“A” and “W” in Table [^respectively. Solid (blue) and dashed 
(red) vertical lines (at —0.52 and —0.13, respectively) are the 
mean values of the distribution functions. Inset: show the 
tails of the distributions near 7 = 1. Error bars (not shown) 
are about the size of the symboles. 


C. The Distribution Function of the End-to-end 
Tangent Vectors 


D. Loop Formation Probability 


In Figure [7j P 1 { r ),L) is plotted against 7 for a DNA 
with L = 36.66 nm(= 99 bp) while its end-to-end distance 
is kept at r = 0. In this Figure, triangles (blue) and 
squares (red) correspond to the models “A” and “W”, 
respectively. As can be seen the two distributions are 
significantly different, but in the vicinity of 7 = 1, they 
asymptotically approach each other, (inset of Figure [7|. 
We found that the peak of the distribution for the AER 
model generally occurs at a smaller 7 compared to the 
WLC model at short length scales (below the persistence 
length). For example, for L = 36.6nm, the most prob¬ 
able values of 7 for the models “A” and “W” are —0.74 
and —0.18, respectively, and with { n f) A = —0.52 and 
(7) W = —0.13 (as indicated in Figure [t] by solid (blue) 
and dashed (red) vertical lines, respectively). This indi¬ 
cates that when the two ends of a stiff chain meet each 
other, the angle between the terminal tangent vectors 
tends to be smaller in the AER model compared to the 
WLC model. This difference reflects the effect of the kink 
formation on the equilibrium structure of the DNA loop. 
The same structure also has been reported by other stud¬ 
ies which are considered the kink in the model 133133. 


Figure [H] compares the the J-factor of the DNA in the 
AER (triangles, blue) and the WLC (squares, red) mod¬ 
els, as obtained in the MC simulations, where filled and 
open sy mbol s co rresp ond to J and Jo, respectively (see 
sections IIB and IIC and equations (111 and (fl2])). The 
solid black curves, are the theoretical predictions for J 
and Jo in the WLC model j3j which perfectly match the 
simulation data. In the case of the AER model, the 
dashed curves are shown as eye-guides. As can be seen, 
at short lengths (below 100bp), the J-factor in the AER 
model (with or with out axial alignment) is several or¬ 
ders of magnitude larger compared to the WLC model. 
As expected, the difference between the two models de¬ 
creases as the DNA length increases, and for length larger 
than the DNA persistence length (~ 150 bp) the models 
are essentially indistinguishable. The same result can 
be obtained by other kinkable models mm- In this 
study we showed that the asymmetry in DNA structure 
may promote the kink formation, in particular, largely 
increases the J-factor at short length scales. As it was 
discussed, the torsional constrain is not considered for 
the both looping probabilities J and J 0 , which leads to 
oscillations of the J-factor as a function of DNA length 
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with a period equal to the DNA helical pitch [3]. 



FIG. 8: (Color online) The J-factor as a function of DNA 
length. MC result for the model “A” and “W” are represented 
by triangles (blue) and squares (red), respectively, where filled 
and open symboles correspond to the cyclization with (J) and 
without (Jo) axial alignment. Shimada-Yamakawa’s theoret¬ 
ical predictions for the J-factor in the model “W” have been 
aslo shown by solid (black) curves [S] - Dashed (blue) curves 
show the trend of the simulation data and do not correspond 
to a theoretical model. Error bars (not shown) are about the 
size of the symboles. 

Recent experimental data of the J-factor of DNA 
molecules, which was performed by Vafabakhslr and Ha, 
have shown short DNA molecules are much more cycliz- 
able than the prediction of the WLC model 0. In this 
experiment the DNA probe was a duplex with two com¬ 
plementary single-stranded overhangs on both ends (two 
sticky ends). Because the single-stranded overhangs are 
10 nucleotide, they are considered as long sticky ends. 
It is expected that joining of such long sticky ends dose 
not require the axial alignment of the duplex ends (Till . 
and the effect of the torsional alignment could be con¬ 
sidered as an oscillation factor in the J-factor. Also they 
can join each other when the end-to-end distance of the 
duplex is less than capture radius, r 0 , which is 10 bp in 
this experiment. We thus evaluate the J-factor with free 
boundary condition and 3.4 nm capture radius for the 
both parameter models, i.e. the models “A” and “W”. 
As Figure [9] shows, while the experimental data signif¬ 
icantly deviate from prediction of the WLC model for 
short DNA molecules, they show a considerable agree¬ 
ment with the AER model at all length scales. The os¬ 
cillations in experimental data is believed to result from 
the torsional alignment between the DNA ends. In this 
issue, it has been suggested that underlying mechanism 
in the case of surface tethered may increase the rate of 
cyclization mm, but this effect can not explain the 
anomalous behaviour of DNA. Other kinkable models 


show a sharp deviation from the WLC model at a critical 
length [42], but the AER model shows a smooth devia¬ 
tion (see Figures [8] and [9]). The length dependence of 
the J-factor in the length range of 57 to 96 bp is much 
weaker in the AER model compared to the WLC model. 
The inset of Figure [9] shows a zoomed view in the length 
range of 50 — 100 bp in full logarithmic plot. To quan¬ 
tify the length dependence we fit power law functions to 
the simulations data points in this range and we found 
Jq ~ L 3 9 and ~ L 16 , where the superscripts indi¬ 
cate the model parameters “A” and “W”. 



FIG. 9: (Color online) J-factor comparison. The triangles 
(blue) and squares (red) correspond to the simulation data 
of the model “A” and “W”, respectively, for arbitrary angel 
between the duplex ends and the capture radius ro = 3.4 nm 
(= 10bp). Circels show the experimental data of Vafabakhsh 
and Ha for surface tethered (filled) and vesicel encapsidated 
(open) experiments [5]. For consistancy with our simulation 
data, we shifed the orginal data by 10 bp to the left for the 
sticky ends. The solid (black) curves are the theoretical pre¬ 
dictions of the WLC model for ro = 0 and 3.4 nm m- Inset 
represents the length range of 50 — 100 bp in full logaritmic 
plot. Dashed (blue) curves show the trend of the simulation 
data and do not correspond to a theoretical model. Error bars 
(not shown) are about the size of the symboles. 


IV. CONCLUSION 

In this paper, we proposed that the asymmetric struc¬ 
ture of DNA can significantly affect the elasticity of DNA 
at short length scales. We have showed that, the extreme 
bendability of DNA at short lengths as well as the kink 
formation in double stranded DNA can originally form 
the asymmetric structure of DNA double helix. To ac¬ 
count for the bending asymmetry we exploited the asym¬ 
metric elastic rod (AER) model, which has been intro¬ 
duced and parametrized in a previous study m- By 









evaluating the effective bending energy and the distribu¬ 
tion function of the end-to-end distance we show that al¬ 
though the AER model is equivalent to a WLC model at 
large length scales, for tightly bent short DNA molecules 
the DNA is much more flexible in the AER model than in 
the WLC model. Using the umbrella sampling method, 
we evaluated the loop formation probability, i.e. the J- 
factor, as a function of the DNA length. We found that 
the unconstrained J-factor in the AER model with cap¬ 
ture radius about 3.4 nm is in excellent agreement with 
the measured experimental data presented in [S, at all 
length scales. This implies that the axial alignment of 
the two ends is not required to join the two juxtaposed 
DNA ends in this experiment. Enforcing an axial align¬ 
ment can induce an 1000-fold change in the J-factor (see 


Figure[7|. This may explain the large dispersity in the ex¬ 
perimental data where DNA molecules with short sticky 
ends are cyclized gnsjcES]- The results presented in this 
paper, suggest that the asymmetric elastic rod model, 
as parametrized in m , is a realistic model to explain 
the elastic behavior of DNA double helix at short length 
scales. 
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